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Chapter 1 


INTRODUCTION 


This report applies a finite volume approach to the two-dimensional 
compressible boundary layer equations. Unlike conventional finite difference 
methods using similarity transformations methods we introduce no equation 
or coordinate transformations of any kind. Instead, we discretize the inte- 
gral boundary layer equations using untransformed variables. We construct 
a suitable computational grid and calculate an initial solution by explicitly 
applying similarity concepts. The resulting numerical scheme is efficient, ac- 
curate, conceptually straightforward, and very simple. Results are given for 
similar and non-similar flows and compared with finite difference solutions 
and data. 


1.1 OVERVIEW 

This chapter presents a background discussion and motivation of the fi- 
nite volume method and a brief summary of the method. Chapter 2 presents 
the governing equations in integral form, non-dimensionalization, and intro- 
duction of a stream function. It is followed by descriptions of the discretiza- 
tion of the governing equations and of grid generation in Chapter 3. A com- 
parison of finite volume and conventional approaches is given in Chapter 4. 
Chapter 5 reviews linearization of the discretized equations and solution of 
the resulting matrix equation. Results for similar and non-similar test cases 
are presented in Chapter 6, and conclusions and recomendations in Chapter 
7. Details of the implementation are given in the appendices. 


8 



1.2 BACKGROUND 


Numerical solution techniques for the boundary layer equations can be 
classified into two categories: integral methods and finite difference methods. 
Integral methods are based on analytic integration of the momentum equa- 
tion (and, for compressible cases, the energy equation) across the boundary 
layer. For example, many integral methods start with the momentum inte- 
gral equation for 2 — Z> incompressible flow without mass transfer: 




( 1 . 1 ) 


where 0 is the momentum thickness, C/ is the skin friction coefficient, and 
H is the ratio of displacement and momentum thicknesses. Equation 1.1 
satisfies the governing equations, not at each point in the boundary layer, 
but in an averaged “global” sense. Since, for a given velocity distribution, 
there are three unknowns in that equation, two additional relations between 
the unknowns must be specified. Equation 1.1 is an ordinary differential 
equation (ODE) in x and is easily integrated numerically given suitable ini- 
tial conditions. Accuracy of the solution depends on the validity of relations 
chosen. 

Finite difference methods model some form of the governing partial dif- 
ferential equations (PDE’s) by difference approximations on a finite grid. 
Because the discrete approximations reduce to the governing equations at 
each grid point as the grid spacing goes to zero, calculations can be per- 
formed to any degree of accuracy. With given initial and boundary con- 
ditions, no further relations or a priori knowledge is required for laminar 
flows. Since initial conditions usually are not available, the usual recourse 
is to transform the PDE’s to so-called similarity form. A typical approach 
for incompressible flow is to use the Falkner-Skan variables 


to give a third order PDE: 


ri = 


frtm + Oifftft) + ^(1 - /^) = <?* — 




y\/Re 
Lg{x) 

\f^-f ^ 

[J'l 


( 1 . 2 ) 


(1.3) 


$ and fj are the non-dimensionalized tangential and transformed normal 
coordinates, respectively. / is the transformed stream function, g(x) is a 
scale factor, and a and are functions of the pressure gradient and g. 
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For similar flow the right hand side of eqn. 1.3 is zero, giving an ordinary 
differential equation (ODE) in the cross-stream direction that is easily solved 
numerically. Similarity form is valid for wedge type flows only, but is widely 
used for generating initial conditions for non-similar flow calculations. 

The primary advantage of integral methods is their simplicity and com- 
putational efficiency. However, the additional closure relations often involve 
empiricism. Also, it is difficult to extend integral methods to problems 
with added physical processes and/or diverse boundary conditions. As a re- 
sult, finite difference methods have become increasingly popular. They are 
extremely accurate for laminar flows and allow straightforward implemen- 
tation of a variety of boundary conditions. However, the transformations to 
similarity variables destroy the physical transparency of the equations and 
are of little utility once the solution has been started. The transformations 
are involved, and the transformed equations are complicated to the point 
that individual terms are not easily identified with the physical process they 
represent. Furthermore, the discrete equations usually can not be shown to 
exactly conserve mass, momentum, or energy, since the governing equation 
is no longer in divergence form. We consider a finite volume approach, be- 
cause that approach has proven to be accurate, yet simple, when applied to 
the Euler and Navier-Stokes equations. 

1.2.1 Finite Volume Method 

Finite volume methods approximate integral forms of conservation equa- 
tions on finite cells (volumes in 3 — D). The equations are approximated by 
summing the fluxes of mass, momentum, and energy from neighboring cells 
into each cell, thereby satisfying the conservation equations over that cell. 
Conceptually, the approach is a blend of the integral and finite difference 
methods. It reduces to an integral method^ if the number of cells is one, 
but is a finite difference approach in the sense that with many cells, the 
governing equations are satisfied locally as well as globally. The method 
has the accuracy and versatility of finite difference schemes, but retains the 
physical transparency of the governing equations. 

A key part of the finite volume approach is the decoupling of the dis- 
cretization of the equations from the grid generation. The equations are 
first discretized for an arbitrary quadrilateral cell in the physical coordinate 
system. The grid must be constructed to provide resolution where needed. 
This is in contrast to transforming the equations to a similarity type coor- 

that method assumes a linear velocity profile. 
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dinate system before discretizing, whereby a rectangular grid is “built into” 
the transformed equations. If a different scaled coordinate system is desired, 
one must introduce a new transformation and discretize the resulting equa- 
tions. With the finite volume approach one needs only to construct a new 
grid and the governing equations are always in a very simple form. 

1.3 SUMMARY 

We present a novel finite volume solution of the two-dimensional com- 
pressible boundary layer equations. A untransformed stream function is used 
to simplify the governing equations, but none of the conventional similarity 
type transformations or scalings are applied. Instead, we use similarity prin- 
ciples explicitly to calculate the initial conditions. The grid is constructed to 
automatically follow the boundary layer growth for both similar and non- 
similar flows. Finite volume discretization of the integral equations gives 
a set of nonlinear equations, which are linearized and solved via Newton’s 
method. Results for several similar and non-similar test cases demonstrate 
good agreement with tabulated solutions, computations obtained with con- 
ventional methods, and experiment. 

The method is applicable to laminar and turbulent flow. We use Suther- 
land’s formula [Schlichting 68] and the Cebeci-Smith algebraic two layer 
model [Cebeci 74] to calculate the laminar and turbulent viscosities, respec- 
tively. 

A variety of boundary conditions can be applied. The method makes no 
particular distinction between solving direct cases where the edge velocity 
is specified and solving inverse cases where the displacement thickness, the 
maiss defect, or some other quantity is specified. 

The novelty of our approach is the use of the integral form of the equa- 
tions without transformation. Instead of incorporating similarity features 
into the governing equations via laborious transformations and scalings, we 
use similarity principles explicitly. The explicit use of similarity is much 
simpler than the use of conventional transformations, but gives equivalent 
results. The finite volume method retains the conceptual simplicity of the in- 
tegral conservation equations without loosing the efficiency and advantages 
of conventional schemes. We recommend it both as a calculation device and 
as a pedagogic tool in presenting numerical boundary layer solutions. 
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Chapter 2 

GOVERNING 

EQUATIONS 


We present the integral form of the boundary layer equations, the equa- 
tion of state, and sample boundary conditions. For convenience we non- 
dimensionalize the variables, but we do not change the form of the conser- 
vation equations. 

The discrete equations in primitive (u, v) variables allow streamwise odd- 
even decoupling of the v velocity profile. We eliminate these oscillations by 
replacing the continuity equation by a stream function, and also use the 
stream function to simplify convective terms in the conservation of mass 
and energy equations. 


2.1 INTEGRAL FORM 

The two-dimensional steady boundary layer equations in integral form 
are: 


^ pu dy — ^ pv dx = 0. 

^ u(pti) dy — ^ v(pu) dx = — ^ P dy — ^ r dx 
u{pH) dy — ^ ^{pH) dx — — ^ q dx 


( 2 . 1 ) 

( 2 . 2 ) 

( 2 . 3 ) 


where x is arclength along and y is the distance normal to a given body. 
Integration is taken to be counter-clockwise, u and v are tangential and 
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normal velocities, P is static pressure, p is density, and H is total enthalpy. 
T and q are the shear stress and enthalpy flux: 


r = (« + (-.) 


(2.4) 


9 = 


.Pri Prt. 


dy 


+ m 




(2.5) 


Subscripts (/) and (t) denote laminar and turbulent quantities, and Pr is the 
Prandtl number ( ~ 0.7 for laminar, and 1.0 for turbulent flows). Reynolds 
averaging ha.s been done for turbulent flow. 

Equations 2.1 - 2.5 are five first order relations for the variables u, v, £T, r, 
and q, /i, />, and P must be given or defined in terms of these variables. The 
laminar viscosity is obtained from Sutherland’s law, and the Cebeci-Smith 
model is used to calculate the turbulent viscosity^. The equation of state 
for a perfect gas gives the density-enthalpy relation. 


7 P 


7 — 1 A * 


( 2 . 6 ) 


Boundary conditions may be applied in the so-called direct or inverse 
mode. In either case, the normal and tangential velocities at the wall must be 
specified; typically, «(a:,0) = 0, and v(x,0) = v«,(a:).* Also, either enthalpy or 
enthalpy Sux at the w'all must be given, as well as the total enthalpy at the 
edge of the boundary layer. “Direct” or “inverse” refers to the specification 
of the remaining (fifth) boundary condition. Direct mode implies specifying 
the edge velocity Uf and amounts to imposing a known pressure field: 


J p 2 

on the boundary layer. Alternatively, one of a number of inverse boundary 
conditions, for example, specified displacement thickness or wall shear, is 
easily applied. 

Initial conditions must be specified for «, v, and H to complete the 
formulation. These will be discussed later. 

Equations 2.1 - 2.5 are a nonlinear parabolic system. They are applicable 
to steady compressible, turbulent boundary layers subject to the assump- 
tions that: 

‘See Appendices B and C, respectively 


(2.7) 
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• normal velocity is small compared to the streamwise velocity, 

• streamwise gradients are small compared to normal gradients, and 

• boundary layer thickness 6 is small compared to body curvature. 

The conservation of mass and momentum equations are in standard form. 
The conservation of energy equation is in total enthalpy (per unit mass) 
form. Its simplicity and its similarity with the form of the momentum 
equation make the enthalpy form useful. Definition of the shear stress and 
enthalpy flux as separate equations gives a first order system of equations. 

We solve equations 2.2 - 2.5 in their most primitive form. After non- 
dimensionalization, a stream function is introduced. No further equation or 
coordinate transformations of any kind are applied. Mass, momentum and 
energy are rigorously conserved in the discrete approximation. 


2.2 NON-DIMENSIONALIZATION 


We non-dimensionalize all variables with the freestream quantities Uqo, 
Poo, and some characteristic length L. The primed variables, below, indicate 
non-dimensional values. 


u' = 


Uo 


v' = 


Uo 


p! ^ 

p 




PooUl 

UooL 

H' = 

H 


q 

OO 

q — 

P.0UI0 

r' — 

r 




PooUio 

fl — 

PooUooL 


X 


y 

JU — 

L 

y — 

L 


( 2 . 8 ) 


With the present non-dimensionalization of /i the Reynolds number does 
not appear as an explicit parameter - it is implicitly contained in the equa- 
tions. The forms of the non-dimensional governing equations, equation of 
state, and viscosity and turbulence models are identical to that of the di- 
mensioned equations. For convenience we drop the primes and henceforth 
refer to the non-dimensional form only. 
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2.3 STREAM FUNCTION FORMULATION 


Equations 2.1 - 2.5 can be discretized immediately. We refer to that 
formulation m the primitive variable formulation since the dependent vari- 
ables include u and v. A disadvantage of the primitive variable formulation 
is its susceptibility to stream-wise odd-even decoupling of the normal ve- 
locity. The disadvantage is purely cosmetic since the cell centered results 
are not affected. However, the odd-even decoupling can be eliminated by 
introducing a stream function into the momentum and energy equations. 
Furthermore, introduction of the stream function into the momentum and 
energy equations zdlows reduction of the number of terms in each equation 
by one, yielding an increase in simplicity and computational efficiency. 

2.3.1 Elimination of the continuity equation 

The definition of the total derivative gives a general relation for ^ = 

p.9, 


If now we define. 



dxl> 

dx 


= -pv 


we recover the standard definition of the stream function 


( 2 . 10 ) 


dip = pudy - pvdx. 

Plugging this relation into the continuity equation, it is evident that con- 
servation of mass is identically satisfied. 

Since normad gridlines are perpendicular to the surface at each stream- 
wise station, dx = 0 along normal gridlines; and, we can replace the conti- 
nuity equation by the truncated definition of the stream function: 


dip = pudyx=eonst.> ( 212 ) 

We apply this definition at each streamwise station, and integrate upwards 
at X = const, to calculate ip. 
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2.3.2 Simplification of momentum and energy equations 

In addition to replacing the continuity equation, the stream function 
can simplify the convective terms in the momentum and energy equations. 
These terms represent the transport of x-momentum pu and total enthalpy 
pH by the u and v velocities. A simple rearrangement of the convective 
velocity and density terms in equations 2.2 and 2.3 gives: 

u • pu dy — V • pu dx u • pu dy — u • pv dx (o 

u • pH dy — V • pH dx ^ H • pu dy — H • pv dx ' ’ ^ 

but, using 2.11, this is simply 


u • pu dy — u ^ pv dx = udtp 
H • pu dy — H • pv dx = Hdxp, 


(2.14) 


2.4 GOVERNING EQUATIONS REVISITED 

For convenience we collect our governing equation set; these are the 
equations that we discretize: 


dtp — P^dyx— const. 


/ ^d^+jpdy + ^ rdx - 0 
j Hdip ^ Q dx — 0 


T = {Pl + Pt} — 

oy 


^ Mi 


[Pri Prti dy 


dH 

-^ + Mi 


r 1 1 5u 
1 — u— 
L Pri\ dy 


(2.15) 


The continuity equation has been replaced by the definition of the stream 
function, and the number of terms in the conservation of mass and conser- 
vation of energy equations is reduced by one. Density is present only in the 
stream function. 
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Chapter 3 


DISCRETIZATION AND 
GRID GENERATION 

I 

[ We begin by discretizing equations 2.15 in finite volume form. This 

form is central to the simplicity of the scheme; conservation is expressed by 
a summation of fluxes around each cell, and each term in the discretized 
equations has a clear physical interpretation. 

Subsequently, we illustrate our method for generating grids for arbitrary 
boundary layers. Streamwise grid lines follow lines of r; = const., where 17 = 
yjh . and cross-stream grid lines are along z ~ const. The parameter h — 
/i(x) is a scaling factor proportional to the boundary layer thickness. With 
this definition of the streamwise grid lines, the computational grid always 
extends above the boundary layer, and grid generation is easily implemented. 

j The final section of this chapter discusses explicit application of similar- 

ity principles to calculate an initial solution. In similar flow each variable 
varies as a known power of x along rf = const. At the initial station, no 
previous solution is available and the number of unknonws is greater then 
the number of equations. We start the solution by replacing the variables at 

t y»+i,y their functional equivalents at y**,;, thus closing the system. Since 

the grid is constructed with streamwise faces along 17 = const. ^ calculating 
a starting solution is very straightforward. 

I 

3.1 FINITE VOLUME DISCRETIZATION 

! 

We adopt a box-type finite volume discretization. All variables are de- 
fined on the four nodes of a control volume (see Figure 3.1). The momentum 
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and energy fluxes through a given cell side are calculated by averaging the 
state variables at adjacent nodes. This flux calculation is second order ac- 
curate on all grids. 

The governing equations discretized in this finite volume form are ex- 
tremely simple. Below are the final forms of the governing equations and 
the corresponding discrete approximations using a box-type difference sten- 
cil. Subscripts i and j refer to node points and numerical subscripts refer 
to cell sides in Figure 3.1. Each cell consists of streamwise faces (sides 1 & 
3) and cross-stream faces (sides 2 & 4). On the latter, Ax is zero and we 
drop the corresponding terms from the equations. 



X 

Figure 3.1: Computational Cell 


Definition of Stream Function: 

dip = pudyj,=con$t. 

Aip2 — P2U2Ay2 = 0 

Conservation of Momentum: 

/ udip + / Pdy + / rdx = 0 

tllAipi + U2Aip2 + UzAips + UiAlpi 

+PiAyi + P2Ay2 + PjAya + PiAy4 

+riAxi + rsAxa = 0 
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Conservation of Energy: 

/ Hdr!>+ § qdx = 0 

+ H2^^2 + -ffsA^s + 

+giAxi + q^Ax^ = 0 

Definition of Shear Stress: 

r = im + Iit)du/dy 
1'2Ay2 - {iH + IH)2AU2 = 0 

Definition of Enthalpy Flux; 

q = [fti/Pr + nt/Prt] dH/dy + m^ [1 - 1/Prj] udu/dy 

q2Ay2 - [lii/Pn + fit/PrtU APj 

—m [1 — l/Pr] U 2 AU 2 = 0 


(3.1) 


In these equations, all ‘‘A ‘‘ terms are simple centered differences in a 
counter-clockwise sense; e.g. 


Ayi = - y<,y 

Ay2 = y»+i,i+i - y»+i,y • 

Cell side values of the functions are linear averages; e.g. 

ui = 0.5 • (ui+ij + Uij) 

U2 = 0.5 • (tXf+lj+l + . 


(3.2) 


(3.3) 


This discretization provides second order accurate approximations on an 
arbitrary grid. It preserves the conceptual simplicity of the governing equa- 
tions, since each term in the discrete equations clearly is the flux of mass, 
momentum, or energy through a cell face. Since it is a Crank-Nicolson 
discretization, the solution is unconditionally stable for linear problems 
[Smith 78]. With all quantities given along line x*, the above equations 
may be solved for quantities along x*+i as discussed later. 


Boundary Conditions 

Boundary conditions are implemented very easily as part of the Newton 
solution procedure. Details of the inclusion of the boundary condition in the 
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Newton matrix will be discussed later. Note here, however, that the discrete 
equations above allow odd-even decoupling of the dependent variables, as do 
most central difference schemes for first derivatives. Thus, it is important 
that the boundary conditions be smoothly varying. To ensure some degree of 
smoothness, we assume boundary conditions are given at cell face midpoints 
and obtain the nodal values by averaging^. 

3.2 SIMILARITY AND GRID GENERATION 

An effective boundary layer grid must follow the boundary layer growth. 
The outer edge of our grids is defined as j/c(x) = f)oo*h{x). h is proportional 
to the boundary layer thickness 5, and we define the constant rjoo so that 
y« > S, Subsequent grid lines follow yj = r)j*hi, where rjj = + and 

the Arjj are chosen to give any desired spacing between the grid lines. We 
choose h{x) to be the incompressible momentum thickness O^y since in most 
laminar and turbulent flows that quantity is proportional to the boundary 
layer thickness [Drela 83]. For similar flows, is known from similarity 
principles. For non-similar flows, we calculate it as a part of the solution at 
each streamwise station. 

3.2.1 Grid generation examples 

The use of similarity principles for grid generation is best illustrated 
by examples. We begin by discussing grid requirements for flat plate and 
stagnation flows, and discuss grid generation for general similar and non- 
similar cases immediately afterwards. 

Grid generation for flat plate and stagnation flow 

The boundary layer thickness on a flat plate with sharp leading edge and 
constant edge velocity is zero at the leading edge and grows as The 

grid must follow this variation. Thus, we define and generate 

the grid below. Because the equations are singular at x = 0, the grid begins 
arbitrarily close to x = 0, but not at x = 0. 

A stagnation flow grid is particularly simple: the boundary layer thick- 
ness at a stagnation point is constant, giving a Cartesian grid. 

^Note: This averaging does not in general preserve similarity. However, since real 
problems are nonsimilar, this is not a shortcoming. Alternatively, a logarithmic average 
could be taken, or, if the data is reasonably smooth, no averaging need be done. 
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Figure 3.2: Computational grid for fiat plate fiow. 


3.2.2 Grid generation for similar and non-similar flows 

We discuss first the problem of generating a grid for any similar flow. 
Here the streamwise variation of 9k is known a priori and grid generation 
is trivial. Subsequently, we discuss the process of calculating $k as part of 
the solution at each station, eis a means of defining the grid in non-similar 
flows. 

SIMILAR. FLOWS 

Grids generated for flat plate and stagnation flows are special cases of grids 
generated for similar flows. An incompressible flow is similar if the edge 
velocity varies eis a power of x: 17* ~ x”*. This is the case in incompressible 
flows around a wedge. In similar flows, the variations of u,r,ip,H, and q 
along r} = const, are, by definition, functions of x only. All of these variables 
and S follow a power law variation in x. Along rj = const, we can write:, 

u = x“«, V = x“^, r = x“% H = x“», q = x“«, and 8 = x“ (3.4) 


A similarity transformation of the differential forms of equations 2.1 - 
2.5 readily yields the desired a's*. If ^ • x is the wedge angle (0 < /? <2),a 
is related to via 


a = 


1 - jg 

2 - ^ ’ 


(3.5) 


^Stagnation and flat plate flow have a«(= m) = 1,0, respectively. 
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and 


o^u = = 1 2a 

ar = aq = 1 — 3a 
a^ = 1 - a 


(3.6) 


For similar flows, we are always given the edge velocity Ue ~ We can 

determine a from ?7g(x), define ~ 0fc;(x,+i/xi)“, and generate our grid 
correspondingly. 


NON-SIMILAR FLOWS 

Method 1. Boundary layer flow is usually not similar. However, if stream- 
wise changes are small, local regions of the flow may be locally similar. 
Method 1 assumes the local flow to be similar and estimates Oki+i> beised on 
a calculated 6ki and a local value of a. 





(3.7) 


This formula reduces to the previous treatment for similar flows, and 
gives a good approximation to &k for non-similar flows. However, if the 
boundary layer thickness grows more rapidly than the quasi-similar assump- 
tion predicts, we must allow the initial edge of the grid, y^, to be far outside 
of the boundary layer, so that the boundary layer always remains within the 
computational domain. 

Method 2. Instead of extrapolating to find this approach calcu- 

lates at each station as a part of the solution. The grid ordinates yi+i,y 
are defined as and is treated as an additional 

unknown in (2.15). To define $k and close equations 2.1 - 2.5 we must add an 
equation - the definition of 0^ - to equations 2.15. This process is discussed 
in more detail in Section 5.1.3. 


3.3 STARTING THE SOLUTION 

At any streamwise station, there are J — 1 discrete control volumes, 2 • J 
nodes (at t t + 1), and five unknowns per node, giving a total of 10 • J 
unknowns. We want to solve 3.1 for the unknowns at station i + 1. When 
the solution at t is known, the number of unknowns decreases to 5 • J. We 
can solve this system since, at each station, we have 5 • ( J — 1) equations 
and 5 boundary conditions. 
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At t = 0, no previous solution is available, and we have only half the 
number of equations needed to solve for the unknowns. Somehow, we must 
generate a starting solution. Our approach is straightforward: 

1. Assume the flow is similar between i = 0 and t = 1 and construct the 
grid accordingly. 

2. Since, in similar flow, each variable V varies like 

along constant rj , replace the unknowns at y»+i,j by those at j/ij with 
the appropriate x-scaling, thereby reducing the number of unknowns 
to 5 • J. 

3. Calculate the terms in the discretized equation, as before. 

Below is an example of the approach applied to incompressible flat plate 
flow. In flat plate flow a„ = 0, a, = -1/2, and = 1/2. To start the 
solution we set at every riji 

Substitutions Fluxes 

«t=i = «,=o => ui= (ti,=o + «.=i)/2 = ti,=o 

V’»=i = V’»=o • ^ V’l = V’»=o • (1 + (a:o/®i)~^^^)/2 

n=i = r,=o • (xo/xi)^/* =» n = U=o’{l + {xo/xiyf^)/2 

etc, => etc, 

(3.8) 

Fluxes through faces 2 and 3 can be defined with similar relations. Now, the 
discretized equations are written in terms of 3 • J variables at i = 0, only, 
and we can solve for the unknowns. The approach is valid for any similar 
flow if the exponent of (xo/*i) is replaced by the appropriate o from (3.6). 

This substitution is equivalent to using our exact knowledge of d/dx 
(fj = const.) to eliminate all streamwise derivatives and to solve the ordinary 
differential equation in q. The key to the procedure is our grid. Since the 
grid lines are generated to follow lines of constant q, the above substitution is 
justified. l>ansformation of the equations to similarity form is not necessary. 

3.3.1 “True” Similarity Solutions 

The use of similarity allows easy calculation of a starting solution. How- 
ever, that solution is not, strictly speaking, a similarity solution. A true 
similarity solution, when scaled by q, must be independent of x. The above 
procedure uses a linear approximation to the streamwise variation of the 
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variables. This is an exact representation of the streamwise terms in stag- 
nation point flows, where all variables are constant or vary linearly in x 
along rf = const. However, when the exponent a ^ 1, 0 the streamwise 
variation is nonlinear and results in a second order truncation error. Thus 
the solution is dependent on x. The error may be eliminated by analytically 
integrating the streamwise terms (see Appendix A) since their streamwise 
variation is known in similar flows. This gives a true similarity solution. 
This correction is of primarily academic interest, however, since the error is 
very small and flows of interest are necessarily non-similar. 
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Chapter 4 


ANALYSIS AND 
COMPARISON 


This chapter compares and contrasts two aspects of the finite volume 
and conventional approaches: 

• explicit versus implicit application of similarity concepts, and 

• conservative versus non-conservative formulations. 

We present a derivation of the Falkner-Skan equation, a popular conventional 
approach, to focus our discussion. 


4.1 CONVENTIONAL APPROACH 


The conventional approach* applies the divergence theorem to the in- 
compressible version of equations 2.15 to get the differential form of the 
governing equations, in so-called divergence form: 


du _ 
dx dy 


(4.1) 


The non-divergence, or non-conservative, form 
du du ,, dUe d^u 


dx dy 


dx 




(4.2) 


(4.3) 


‘The derivation in this section closely follows [Schlichting 68] (Ch. 8 & 9). 
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is obtained by differentiating the products in the momentum equation, in- 
troducing the Bernoulli equation, and simplifying. The continuity equation 
can be integrated by introducing a streamfunction: 


u 


dip 

dy 


V = — 


dip 


Introducing (4.4) into equation 4.3 gives 
dip d^ip dip d^ ip 


dy dxdy 


= 

dx dy^ ® dx dy^ 


(4.4) 


(4.5) 


The continuity equation hzis been eliminated at the cost of raising the order 
of the governing PDE by one. 

Equation 4.5 can be cast in similarity form by application of one of a 
number of possible similarity transformations^. The Falkner-Skan transfor- 
mation employs independent variables ^ and rj 


X 

L ’ 

and a dimensionless stream function 



(4.6) 


fU>v) = V’(®,y) 


s/r 

LU,{x)g{x) • 


(4.7) 


ii is a Reynolds number, and L is a reference length. The scale factor flr(x) 
will be determined later. 

Introduction of (4.6) and (4.7) into the momentum equation 4.5 yields, 
after some algebra, the governing third order PDE for boundary layer flow: 


fqrtrj + frjrt + /n) “ jj 9^ 








nn 


dj; 

ae. 


where ol and p are: 




(4.8) 


(4.9) 


Equation 4.8 is the counterpart of the incompressible form of equations 2.15, 
and is the starting point of many conventional numerical analyses^. 

^We show later that a valid similarity transformation for the boundary layer equations 
must reduce to one general form. 

®Most numerical approaches now take the somewhat circuitous route to convert (4.8) 
back into a system of three first order PDE’s in { and ij before discretization and solution. 
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Similar solutions to equation 4,8 exist if /, /rj, B.nd a and (3 do not depend 
on z. Then the right hand side of (4.8) vanishes and an ODE in r) remains 
on the left hand side. The constraint on a and gives two equations for the 
edge velocity and the scale factor. After some algebra one obtains 


9 = 


2 X Uoo 
m + lL Ue ’ 


(4.10) 


and 

tTe = C-i"‘. (4.11) 

Potential velocity distributions with this form are found in wedge flows. 


4.2 EFFECT OF SIMILARITY 

If Ue = Cx”*, the flow is similar and equation 4.8 reduces to an ODE, 
allowing straightforward calculation of a starting solution. Although the 
flnite volume approach uses no transformation, similarity principles can be 
used to stMt the solution. The two sections below discuss the parallels 
between the two approaches. For the purposes of this analysis, constant 
factors may be ignored. 


4.2.1 Conventional Approach 

For similar flows the scale factor g is 

_ j 2 X 
^ y m + lL Cx”» 

= , / - _ 

Y (m+ l)CZr 

Thus, the independent variable rj becomes 


_ 

(C,i = constant), and the dependent variables are 

V'(®,y) = \j — /(»?) 

V m + 1 


(4.12) 


(4.13) 


(4.14) 
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u(®,y) = Cx”"/'(»j) 


(4.15) 

(4.16) 


v{x,y) = 



m “ 1 9/‘ 
m + 1^ 9 t7. 


Thus, in similar flow the stream function and the velocity components^ are 
known functions (up to a constant multiplicative factor) of x along t] = 
const. 

The "character” of the Falkner-Skan transformation becomes clear if 
we consider the dependent variables 4.14 - 4.16 for constant rj. Then, the 
relations become 

V>(x) =C^x(”‘+^)/2 

u(x) = CuX^ (4.17) 

v(x) = C’„x (”‘-^)/2 

where C^, Cu, and depend on r). The streamwise and normal velocities 
are of no immediate interest since they are eliminated by the stream function. 
Consider the transformed streamfunction, however. From equation 4.7 


= V’(®.y) 


yjR 

LU,{x)g{x) • 


Inserting the definitions for and ^(x) gives 


/(C,»?) = V’(®,y) 


^ 


(4.:s) 


Comparison with equation 4.17 shows that, in similar flow. / = const. 
along constant r} . This certainly is not surprising. It is the crux of the 
similarity transformation. If it were not the case, equation 4.8 would not 
reduce to an ODE in rj - the solution would not be similar. 

Any number of similarity transformations are possible (and a large num- 
ber have appeared in the literature). However, for similar flows they must 
have two points in common: 


1. The ordinate must scale as g{x). 

^Note that the normal velocity in the conventional approach is defined by an ODE in 
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2. The transformed variables along any tj = const, line must be constant. 
Thus, the transformed variables (’) must be of the form (from (4.17)) 




j.(m+l)/2 

U 

V 


(4.20) 


In this limited sense, similarity transformations of the incompressible bound- 
ary layer equations differ only in the degree of complexity that hides (4.17). 

A noteworthy transformation in this regard is that due to Drela [Drela 83] . 
Drela writes the governing equations in first order form and transforms each 
primitive variable individually by the appropriate combination of 17* and 
j(x). For example, ix' = u/Ue and rp' = \p / {U tg{x)) . In similar flow, 
Ug ~ X*", Ugg{x) ~ etc. and his governing equations reduce to 

ode’s, as they must. His approach is one step towards simplifyng the gov- 
erning equations - each transformed equation still represents a conservation 
law. However, the transformations and manipulations are still involved, the 
number of terms per equation is increased, and each term does not have a 
physical meaning. 


4.2.2 Finite Volume Approach 

The finite volume approach leaves the equations intact; rather, the grid 
is scaled to adapt to the flow. The normal location of the grid lines lie along 
y = fjh(x). With h{x) defined as the incompressible momentum thickness, 
h{x) is equivalent to y(x) since in similar flow 9k S ^ (cf. eqn. 

4.12). 

Having defined a computational mesh with grid lines that lie along con- 
stant ry, each of the similarity relations 4.17 holds true. To calculate a 
starting solution, we express the variables at the unknown upstream station 
in terms of the values at the assumed known station (see Sec. 3.3), and ap- 
ply the Newton procedure until convergence. In effect, the discrete PDE’s 
reduce to ordinary difference equations. 

The conventional and finite volume approaches can be thought of m, 
respectively, applying similarity principles explicitly (by transforming the 
equations) and implicitly (by scaling the grid). The grid scaling in the finite 
volume approach is functionally equivalent to the similarity transformation 
used in conventional schemes but is much simpler. 
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4.3 CONSERVATION 


Conservation of mass, momentum, and energy means that in the discrete 
approximation to conservation equations, no spurious mass, momentum, or 
energy is created. Conserving mass, momentum, and energy is probably not 
mandatory in the subsonic region of the boundary layer, but is nevertheless a 
good feature. However, shocks may penetrate the essentially inviscid region 
near y = and strict conservation there may be essential for accurate 
prediction of shock strength and location. Conservation is one feature of 
the finite volume approach; however, conservation can not be shown for the 
discrete approximation to the transformed governing equation 4.8. 

The finite volume method approximates the integral conservation form 
of the boundary layer equations on finite grid cells. The laws of conservation 
of metss, momentum, and energy are approximated by summing the fluxes of 
those quantities into and out of each cell. If one takes a line integral of the 
fluxes around an arbitrary number of cells, or, equivalently, sums the line 
integrals of each of the enclosed cells, the only contribution to that integral 
is due to fluxes through faces at the boundary of our contour. The fluxes on 
the inside of the contours cancel, since the fluxes leaving each cell are just 
those entering the surrounding cells, and every flux is accounted for. 

Differential forms of the governing equations can be shown to be conser- 
vative if they can be written in the divergence form 


dx dy 


0 . 


(4.21) 


and if the equations are discretized consistently. To show this, we sim- 
ply apply the divergence theorem to equation 4.21 and obtain (assuming 
counter-clockwise integration) 


^ Fdy - ^ Gdx = 0 , (4.22) 

and show conservation as outlined above. 

Equations 4.1, 4.2, and 4.21 are in this form; however, equation 4.8 is 
not. Since there seems to be no direct means of writing (4.8) in divergence 
form, we can not show discretizations based on (4.8) to be conservative. 
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Chapter 5 


LINEARIZATION AND 
SOLUTION 


This chapter discusses linearization of difference equations 3.1 and so- 
lution of the resulting matrix of linear equations with Newton’s method. 
Because of the compact nature of the discretization in Section 3.1 the lin- 
earized system of equations forms a block tridiagonal matrix which can be 
inverted very efficiently. Sections 5.2 and 5.3 discuss the setup and inver- 
sion of this block matrix. Section 5.2 also generalizes the right hand side of 
this matrix equation to include global unknowns, which allow us to treat a 
variety of boundary conditions easily. 

Much of the material on linearization, matrix inversion, and treatment 
of arbitrary boundary conditions is discussed in bits and pieces elsewhere. 
We present it here for the sake of completeness. 


5.1 NEWTON LINEARIZATION 


Equations 2.15 and their discretized counterparts 3.1 are a non-linear 
system of equations. Algebraic inversion to find the state variables is im- 
possible. Instead, we apply Newton’s method: we linearize the non-linear 
system, and invert the resulting linear system successively until equations 
3.1 are identically satisfied. We review a one-dimensional application of this 
method briefly, to motivate its extension to multi-dimensions. 
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Figure 5.1: Newton Iteration in 1-Dimension 

5.1.1 Rootfinding in 1-Dimension 

In one dimension, we have a scalar function / = /(v), along with a test 
solution V = v*. In general /(v") ^ 0, therefore we seek a root v* so that 
/(v*) = 0. Newton’s method constructs the truncated, two term Taylor 
series: 

/(v) = 0 = /K) + /V)Av, (5.1) - 

for /(v), about w®. Subsequent approximations are found from = 
v‘ + Au, where, from (5.1), 

= -4^ (5.2) 

Geometrically, this process consists of calculating the tangent of / at 
f{v*)y and sliding down it to (see Fig. 5.1). We repeat the process until 
Av is smaller than some predetermined value. 

If / is well behaved, and we are far from local maxima or minima, each 
iteration gives a much improved V 2 due of v; therefore, we quickly converge 
to the solution. If f{v) is a linear function of v, /' = const, and one ap- 
plication of formula 5.2 immediately yields the desired result. If, however, 
/(v) has local maxima or minima near v, f\v) 0 and the algorithm fails 
spectacularly. 

5.1.2 Rootfinding in iV-Dimensions 

The extension of Newton’s method from 1 to J-dimensions is conceptu- 
ally straightforward. Now we search for the solution vector V for the system 
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of equations F(V) = 0. Instead of a single Taylor series expansion in one 
variable, J expansions in J variables are required for a vector of unknowns 
AV. Instead of the ordinary derivative /', we calculate J x J partial deriva- 
tives, dfi/dvji,j = !••• J, and / can depend directly or indirectly on V. 
This results in the system of J equations, below: 

/i(V) = 0 = /x(V‘) -h fA Avi -h fAAuj 4- |A Avs -h ... 

/2(V) = 0 = /2(V) + §A Avi + §Aav 2 -h fAAua + ... 

/S(V) = 0 = /3(V‘) + §AAVx + ||AV2 + |AAv3 + ... 


fj(V) = 0 = /y(V‘) -h |AAvi -h §Aav 2 -f- ^Av3 + ... + I^Avj 

(5.3) 

The process of iteration is the same: calculate successive AV's and 
iterate to drive F(V) = 0. However, instead of sliding down lines to get 
the V* that gives the intersection of /(v‘) with 0, we now slide along J 
dimensional planes (hyperplanes) in search of their intersection. Also, the 
AV's can no longer be found by simple division. Inversion of the system to 
find AV is a major task and is discussed in Section 5.3. Systems of specific 
interest to us are given below. 


Incompressible laminar boundary layer flow In incompressible flow 
there are three equations in three unknowns at each cell j, giving (for J — 1 
cells) a system of 3(J - 1) equations in ZJ unknowns^. The resulting linear 
system for cell j is 


Rn.,(V) = 0 = Rm,M + 
fi*,(V) = 0=ii*,(V<) + 
R,.(V) =0= Ji,,(V‘) + 


dR, 


dR„ 


dR„ 


anf'^^3 + + ( • 

aRl. . dR^, . , dR^, , , 

-I- s^A^j + -5;^ Ary + ( . 

lUf Any + -g^AV'y + -y?fAry -h ( . 


• )i+i 

• )i+i 

• )i+i 
(5.4) 


Rmp Riij, and I?ry are the residuals of the momentum equation, and the 
definitions of stream function and shear stress, respectively (equations 3.1). 
The unknowns are AV = [Au Afp Ar at j and j + 1. V*' is the vector 
of initial guesses. 


*If the solution at t is known, each residual 3.1 depends on u, ip, and r at (»-|- l,j) and 
(t + l,j+ 1). However, the unknowns are shared by adjacent cells leading to a total of ZJ 
variables. Three boundary conditions must be specified for closure. 
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Compressible turbulent boundary layer flow In compressible flow, 
there are five dependent variables, and we linearize the five discretized equa- 
tions with respect toV=[u^rjffg]^. For compressible and turbulent 
flows the residuals are implicit as well as explicit functions of V. This dual 
dependence of the residuals on V must be accounted for by applying the 
chain rule of diflFerentiation when appropriate. 

For example, the residual of the stream function for compressible de- 
pends on u and x() explicitly, and on u and H implicitly. The implicit de- 
pendence is due to the equation of state: 


^ 7 P 

Application of the chain rule gives (for a given P) 


du 

dR 

dH 


dR dp 
dp du 

dR dp 
Jp'dH 


dR 7 

dp-f-1 



dR 7 -P 


(5.5) 


(5.6) 


Similar cross-differentiation should be performed for all implicit functions of 
ti, r, and q to ensure quadratic convergence of the Newton algorithm. 

The cost of neglecting the variation of any quantity is high. Figure 5.2 
shows the number of iterations needed for convergence when turbulent vis- 
cosity terms are linearized and not linearized. For higher Reynolds numbers, 
each station will require even more iterations. 


5.1.3 Global unknowns 


It is advantageous to include in the set of variables two global unknowns, 
A$ic and AU^. We refer to A0k and AU^ as global unknowns because they 
are unknown quantities that are functions of x only, and do not vary across 
the boundary layer. By allowing A9jc to be an unknown, we calculate the 
grid as part of the solution at each station {yij = usually is 

known a prioriy and AU^ = 0. However, in inverse calculations a quantity 
other than Uc is specified and AUc must be treated as an unknown. 

The global unknowns add two terms to each of equations 5.3, 


A/, 

dOk 


and 


dR 

dU, 


AUe 


(5.7) 


I 
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Figure 5.2; Effect of linearizing eddy viscosity on number of iterations per 
streamwise station. 

Although these terms are added at esu:h cell, the total number of variables 
increeises by only two since A9k and AUg are constant in the cross-stream 
direction. Two corresponding relations must be introduced to close the 
system of equations. The (incompressible) definition of 9 k, 

upon discretization gives one relation. The second must be derived from the 
fifth specified boundary condition. If Ue is specified, 

AUg = Ugpgg - Ug (5.9) 

is sufficient. If some other (inverse) boundary condition is specified, the form 
of (5.9) must be derived from the definition of that boundary condition (see 
Appendix D). 

The turbulence model contains two additional global unknowns, velocity 
defect Au and shear velocity Uf For simplicity, we do not treat them as such. 
In effect, this causes the influence of their variation on the residual to lag by 
one iteration, and the convergence rate may decrease. However, the velocity 
defect typically varies slowly &om station to station; thus, neglecting its 
variation has little effect on the convergence rate. The global changes in the 
shear velocity can be adequately modeled by the local variations in shear 
stress and density since Ur = f{r^,Pw) and ot pj in the inner 

layer. 
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5.1.4 Boundary Conditions 

Applying the Newton method to the set of five equations over J — 1 cells 
gives a set of 5 • ( J — 1) linear equations (5.3) at each streamwise station, for 
5 • J + 2 variables (if the solution at the previous station is known). The two 
auxiliary relations and five boundary conditions supply the equations needed 
for closure. The boundary conditions are implemented by initializing the 
applicable state variable to its boundary value, and requiring all subsequent 
changes to that variable be zero. This gives five relations of the form 

Au?,. = 0 

(the no slip condition in this example) that close the system. 

5.2 BLOCK MATRIX SETUP 

The system of equations 5.3 is complete when boundary conditions, 
global unknowns, and auxiliary relations are included. With two global un- 
knowns, a computational region of J — 1 cells results in a (5 J + 2) X (5 J + 2) 
matrix: 

Ai Cl Gi 

B2 A2 C2 G2 

Bj Aj Gj 

ZX XX XX XX XX 

XX XX XX XX XX 

AVj is the vector of unknowns [Au Atp At AH Aq Ajy Bjy and 
Cj are 5x5 matrices containing boundary condition terms and the deriva- 
tives dRj/ dVj and dRj/dVj^i. Rj is the negative of the vector of residuals of 
the five equations at cell j, while Gj and Hj are vectors containing d'Bj/dOk 
and respectively. The (J + l)’th and (J + 2)’th line each contains 

one scalar equation, the linearized momentum thickness (5.8) and (5,9) (or 
some other specified boundary condition), respectively. Derivation of these 
relations is discussed in Appendix D. 

Due to the compact discretization, each residual depends only on vari- 
ables at the neighboring node points and on the global unknowns. With 
imposition of the boundary conditions, the structure of the 5 J X 5 J sub- 
matrix (sans global terms) is blocktridiagonal. One may take advantage of 




(5.10) 
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this structure by transposing the columns contaifting the global terms to the 
right hand side of (5.11) and treating the global unknowns as quantities to 
be determined after inversion. The resulting matrix of equations is 


■ Ai Cl 
B2 A2 C2 


AVi' 

AVa 


Ri' 

Ra 


Gi‘ 

G2 

-AUe 

Hi' 

Ha 

Bj Aj 


AVj 


.Rj. 




Hj. 

(i 


All matrices and vectors in equation 5.12 are listed in Appendix E. Equa- 
tions J -h 1 and J -f 2 will be reintroduced after this system of equations is 
inverted. 


5.3 BLOCK TRIDIAGONAL INVERSION 


Solution of the block tridiagonal matrix is similar to solution of a scalar 
tridiagonal matrices. The algorithm is the same, with scalar operations 
replaced by matrix operations. At each iteration, we form the recurrence 
relations 

Di = Ai 

Ej = [I>7^]Cy l<j<J (5.13) 

Dj = Aj - 2<j<J. 

For 3x3 blocks we form the Inverse via Cramer’s rule. This is inefficient for 
5x5 blocks, where instead we factorize the matrices by Gauss elimination. 
Applying the same operations to the right hand side gives 


Zi = 


z< = If 7‘IIR,- - 

2 < i < J 

zi, =[fr‘iGi 


Zl,- = - |B(JZ1,-,1 

2 < i < J 

Z2, =Pr'lHi 


Z2, = - |B,lZ2y.J 

2 < i < J 

left and right hand sides gives the 

intermediate result 


(5.14) 


I El 
I E2 



AVi' 

AVa 


Zl' 

Z 2 

-A0k 

'Zli' 

ZI 2 

-AUe 

'Z2i' 

Z2a 


AVj 


.Zj. 


.Zlj. 


_Z2j 


(5.15) 
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where I is the identity matrix. 

Finally, we eliminate the upper band of E matrices: subtract row 

j from row — 1 to get 


■ AVi ■ 
AV 2 

— 

1 

1 

-A0k 

■ ZZli ■ 
ZZ 12 

-AUe 

■ ZZ2i ' 

ZZ22 



.ZZj. 


ZZlj 


ZZ2j 


where, for example, 


ZZj =Zj 

ZZy_x = Zj — Ej^iTjZij 2 < j < ,7, 

etc. 

At each j, we now have the scalar equations 


(5.17) 


Auy = Z Z\j — Adjt • ZZl\j — AC/^g • ZZ2ij 

A^jf — ZZ^j — A^jt * ZZI2,) — A£/^g ‘ ZZ2^j 

Ary = ZZsj ~ A^j. • ZZl^j — AC/^g * ZZ23j (5.18) 

AHj = ZZ^j - A 9 k • ZZUj - AUe • .^.^24.,- 

Aqj = ZZ5J - A 0 k • ZZUj - AU, • ^^25’,,- 

that must be solved for Auy, AV»y, Ary, AHj, and Agy. ZZ, ZZl, and ZZ2 
are known quantities. The two auxilliary relations determine A 0 k and AUg. 
We write them in the form 


AiA0* + BiA17g = i?i rt; 10^ 

AiA0k + B2AU, = R2 , ^ 

(see Appendix D.2) calculate A0k and AUe using Cramer’s rule, and, sub- 
sequently, calculate AVy for all j. 

A convenient convergence criterion is the size of (AV /V)y. If mai( AV /V)y < 
c, the Newton iteration has converged to a solution at that streamwise sta- 
tion. If not, we update = VJ- — AVy and repeat the process. 


5.3.1 CPU Requirements 

Solution of the block tridiagonal Newton matrix leads to very modest 
computer time requirements; especially, if the inversion procedure is op- 
timized by taking zero columns and rows in each submatrix into account, 
For example, solution of the incompressible equations with the Cebeci-Smith 
turbulence model and 40 cells across the boundary layer takes approximately 
0.12 seconds/iteration of CPU time on a Digital Microvax iP. Typically, 2—6 

^Comparable to the Digital VAX 11/750 
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iterations are required for convergence. Thus, a turbulent test case with 50 
stations and 40 cells across the boundary layer requires well under 30 seconds 
of CPU time on the Microvax IL 
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Chapter 6 

RESULTS 


The finite volume scheme (FVS) has been applied to a number of test 
cases with excellent results. We present several typical similar and non- 
similar test cases below. Results obtained for similar fiows are compared to 
solutions obtained with Keller’s Box Scheme (KBS) [Keller 70]. Non-similar 
results are compared to solutions obtained with the method of Drela^ (DS) 
[Drela 83]. Results for turbulent flow cases are compared with experiment, 
and, in one case, with a previous solution. 

6.1 SIMILAR FLOW 

Similar flow, though of little practical interest, provides ideal test cases 
for comparison of the finite volume scheme with conventional methods. The 
flat plate non-dimensional wall shear stress 

( 6 . 1 ) 

and displacement thickness 



values, with the percent error are given in Table 6.1. They have been calcu- 
lated with KBS and FVS for various grid spacings^. The upper boundary 

^Drela’s method is a conventional (by our def.) second order accurate scheme that has 
been validated on a number of test cases. 

^Solutions were “optimized” by aposienori choice of the 6 (6 = Arjj^i/ Arij) that gave 
the most accurate solution for each method and case with a ten cell base grid. Denser 
grids were generated by subdividing each base grid cell into equal intervals. 
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of the computational grid was set at 


for both schemes^. 


Table 6.1: Flat plate wall shear/percent error: 



WALL SHEAR 

DISP. THICK. 

# cells 

Keller 

Fin. Vol. 

Keller 

FIN. Vol. 

b 

1.0 

1.2 

1.0 

1.2 

10 

0.33085 

O.S6S% 

0.33700 

■1.488% 

1.77121 

■2.982% 

1.69312 

1.605% 

20 

0.33177 

0.087% 

0.33328 

-0.868% 

1.73342 

■0.787% 

1.71353 

0.419% 

40 

0.33199 

0.022% 

0.33237 

■0.092% 

1.72393 

■0.185% 

1.71899 

0.102% 

80 

0.33204 

0.006% 

0.33213 

■0.028% 

1.72156 

■0.047% 

1.72037 

0.022% 

Exact 

0.33206 

1.72074 


KBS wall shear results are somewhat more accurate than FVS shear 
results, but the finite volume scheme predicts the displacement thickness 
slightly more accurately than does KBS. Both sets of results are extremely 
accurate, however, and the differences are minute. The errors for both 
schemes decrease as the square of the grid spacing, indicating second order 
accurate discretizations. 

Table 6.2 gives the non-dimensional shear and displacement thickness 
for stagnation point flow as calculated by FVS and KBS. Definitions (6.1), 
(6.2), and (6.3), have been used with Uoq/x replaced by dU^/dx and tjoo = 4. 
The schemes predict the wall shear equally well, but FVS again gives more 
accurate displacement thicknesses. This might be expected from the basic 
formulation. The KBS approximates the differential equation and therefore 
gives greater accuracy for pointwise quantities. The FVS approximates the 
integral equations thus yielding greater accuracy of global quantities. 

®Note: this if differs from the if in if = yjd by a factor of 0.664 
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Table 6.2: Stag, point wall shear/perceni error: 



WALL SHEAR 

DISP. THICK. 

# cells 

Keller 

Fin. Vol. 

Keller 

FIN. Vol. 

b 

1.2 

1.2 

1.2 

1.2 

10 


1.237747 

■ 0 . 418 % 

0.65871 

■ 1 . 666 % 

0.64426 

0.563% 

20 

1.23126 

0.108% 

1.23388 

-0.105% 

0.65059 

■ 0 . 414 % 

0.64701 

0.139% 

40 

1.23226 

0.027% 

1.23291 

■0.026% 

0.64853 

■0.095% 

0.64772 

0.029% 

80 

1.23251 

0.007% 

1.23267 

■0.007% 

0.64801 

■0.015% 

0.64791 

0 . 001 % 

Exact 

1.232588 

0.64791 


6.2 NON-SIMILAR FLOW 

6.2.1 Laminar Flow 

Figure 6.1 shows the distribution of momentum thickness along an infi- 
nite cylinder as calculated by FVS and DS with the potential flow pressure 
distribution. Fifty points in the normal direction and Ax = tt/ 180 were 
used. This test case is a good example of non-similar laminar flow around a 
blunt body. At small 0, u increzises linearly and the flow approximates stag- 
nation flow. After 9 — 90® the velocity decreases and the pressure gradient 
becomes positve; and, finally at ^ = 106® — 109®, the flow separates. 

The two curves are almost identical. In both cases, the solution di- 
verged at 0 = 107 indicating separation of the flow. Figure 6.2 gives the 
corresponding wall shear stresses. Again, the curves are almost identical. 

Howarth’s flow is a popular, if artificial, laminar test C 2 tse. It is defined 
by a linearly decreasing edge velocity: Ug = 1 — x/8. This definition gives a 
singular solution at x = 0, an increasingly positive pressure gradient at larger 
X, and separation at x .96. In this calculation there are 100 cross-stream 
points and Ax = 0.01. Figures 6.3 and 6.4 show the momentum thickness 
and wall shear stress, respectively. The curves are nearly indistinguishable 
away from the singular point. Both solutions fail to converge at x = 0.965, 
indicating separation at about x = 0.96 (the previous station is at 0.955). A 
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Figure 6.2: Wall shear on infinite cylinder 
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Figure 6.3: Momentum thickness in Howarth’s flow 



^.OOOE+00 0.250 0500 0.750 ^00 


X 

Figure 6.4: Wall shear in Howarth’s flow 
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grid density study in [Loyd 85] showed no significant changes as the number 
of grid points was decre 2 ised to 25 x 25. 

6.2.2 Turbulent Flow 

Finally, we compare numerical results of two turbulent flow cases with 
experiment. The primary purpose of these test cases is to demonstrate the 
straightforward extension of the finite volume scheme to turbulent flows. 
As such, no attempt was made to “fine tune” the turbulence model or to 
“smooth” either input data or numerical results to attain better agreement 
with experiment. 

In the following figures unconnected symbols represent either the experi- 
mental data or previous calculations, while connected symbols are the finite 
volume results. All test caises were run with 40 cells across the boundary 
layer and a stretching ratio b — 1.15. 

Flows 2100 and 2400 from the 1968 Stanford Conference [Coles 69] are 
fully turbulent, and essentially incompressible. Initial conditions must be 
assumed to start the computations. A possible approach is to use the initial 
f7— velocity profile data, integrate it to get the appropriate tpij, and input it 
to the turbulence model to calculate the apparent stresses and, subsequently, 
the shear stresses. However, this method is difficult to implement if few 
profile points are available, since interpolation is necessary and a suitable 
near wall model must be used. Instead, our approach hzis been very simple. 
We assume that the fully turbulent flow is “similar” in a small region near 
the first station, and choose the exponent in Uc to get the best possible 

agreement of the “similarity” solution with the data at the first station. 

Flow 2100 is the boundary layer on a large airfoil-like body, at Reynolds 
number^ Re = 5.0 • 10®. The edge velocity distribution is shown in Figure 
6.5. The pressure gradient is negative at first, then nearly constant, and 
finally positive. The experimental data indicates the onset of separation 
at the final stations. Data is available at 38 streamwise stations, but at 
less than ten cross-stream data points near the leading edge. Calculation 
stations were taken at the locations of the experimental data. 

Figure 6.6 shows the measured and calculated mean velocity profiles at 
the first three stations. The profiles were non-dimensionalized by the edge 
velocity at the first station. For clarity, the profiles have been plotted for 
^/Uci > 0.5. The profiles at all three stations agree quite well with the 

^Based on reference quantities at first station. 
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0.00 3.75 7.50 11.3 15.0 18.8 22.5 26.3 30.0 

X(ft) 

Figure 6.5: Flow 2100: Measured edge velocity 

experimental data, thus validating the assumptions made in calculating the 
initial conditions. 

Figure 6.7 gives the experimental and predicted skin friction coefficients 
for Flow 2100. The predicted values are in excellent agreement with the 
data for x < 20ft. For x> 20ft. (approaching separation) the computation 
overpredicts the skin friction. 

The shape factor jff(x) is plotted in Figure 6.8, and the momentum 
thickness in Figure 6.9. Again, the trends and variations of the experiment 
are modeled well except near the separation point. Inability of the turbu- 
lence model to accurately predict the eddy viscosity near separation is the 
probable reason for this discrepancy. 

Flow 2400 begins in a fairly strong positive pressure gradient, which 
abruptly decreases to zero (see Fig. 6.10). Extensive cross-stream data is 
available, but only 7 streamwise stations were measured. Twenty stations 
were used in the computation. The edge velocity at the intermediate stations 
was obtained by linear interpolation of the experimental data and is given 
in Figure 6.10. No attempt was made to smooth out the discontinuity in 
dUjdx. 
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Figure 6.6: Flow 2100: Mean velocity profiles 


UJ 



Figure 6.7: Flow 2100: Computed and experimental skin friction 
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0.00 3.75 7.50 11.3 15.0 18.8 22.5 26.3 30.0 

X(ft) 

Figure 6.9: Flow 2100: Computed and experimental 6 
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3.00 3.75 4.50 5.25 6.00 6.75 7.50 8.25 9.00 

X(ft) 

Figure 6.10: Flow 2400: Measured edge velocity 

Due to the adverse pressure gradient the numerical results are very sen- 
sitive to the assumed power variation of Ue at the first station. Velocity 
profiles at the first three stations are shown in Figure 6.11. The agreement 
is not as good as in the previous case. 

Figure 6.12 gives the predicted and measured skin friction distribution. 
A maximum error of approximately 10% is evident at the first and second 
stations. The influence of the initial conditions decreases as the solution 
marches downstream and the values of the predicted skin friction approach 
the experimental data. • 

Figures 6.13 and 6.14 show the shape factor and the momentum thick- 
ness, respectively. The magnitude of the error in the predicted values par- 
allels that of the skin friction. 

The results for turbulent flows are comparable in accuracy to finite differ- 
ence calculations in [Coles 68] computed with the Cebeci-Smith turbulence 
model. In the present study we have attempted only to show that a finite 
volume method is as accurate as corresponding finite difference methods. 
We have not attempted to improve the turbulence model. 

Finally we present results for a compressible flow case. It is the tur- 
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Figure 6.11: Flow 2400: Mean velocity profiles 



X(ft) 

Figure 6.12: Flow 2400: Computed and experimental skin friction 
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x/c 

Figure 6.15: Circular arc airfoil: Pressure distribution 


bulent boundary layer on a bicircular arc airfoil with 18% thickness ratio, 
at Moo — 0.7425, Rsoo = 4 • 10®, and zero angle of attack. This flow case 
was previously calculated as a viscous/inviscid coupling problem by Wigton 
and Holt [Wigton 81]. They used a potential solver in the inviscid region, 
and Green’s Lag entrainment scheme in the boundary layer. A plot of the 
pressure distribution calculated with their coupling scheme is shown in Fig. 
6.15. The flow accelerates around the leading edge, goes through a weak 
shock at about 70% chord, decelerates, and separates just upstream of the 
trailing edge. 

Specification of the edge velocity (direct mode) in the region of separated 
flow leads to an ill-posed set of equations. This problem can be circumvented 
by applying inverse boundary conditions to the problem for example by 
specifying the wall shear or displacement thickness®. Also, to satisfy zone 
of dependency principles, the FLARE technique was implemented® We pose 
the problem as a mixed direct /inverse boundary condition problem, with 

^This is discussed in detail in Appendix D .2 

®...by setting f ud^2 = / udip4, = 0 for U2 + U4 < 0 . However, the energy equation was 
left unaltered since implementation of the FLARE technique there led to oscillations and 
instability in the total enthalpy. 
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Figure 6.16: Circular arc airfoil: Edge Velocity 


the edge velocity specified up to 68% chord, and the displacement thickness 
specified subsequently. 

Some eissumptions about the transition from laminar to turbulent flow 
must be made since [Wigton 81] does not discuss this issue. However, all 
figures of calculated quantities in [Wigton 81] begin at x/c = 0.2. Thus, 
we arbitrarily assumed the flow transitions to turbulence linearly between 
x/c — 0.1 and x/c — 0.2. 

Figure 6.16 gives the edge velocity. Wigton’s data is shown as uncon- 
nected symbols, while finite volume results are plotted as connected symbols. 
The edge velocity due to Wigton et al. was calculated from the reported 
pressure distribution (Fig. 6.15) by assuming constant stagnation pressure. 

Figure 6.17 shows the displacement thickness as calculated (upstream 
of x/c = 0.68), and the results taken from [Wigton 81]. Downstream of 
x/c = 0.68 the two curves coincide since here displacement thickness was 
specified, and the agreement upstream of that point is excellent. 

Figure 6.18 gives the finite volume skin friction results as well as the skin 
friction given in [Wigton 81] for x/c > 0.2. The agreement with Wigton’s 
results is excellent up to about 70% chord, however the predictions of the 
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0.00 0.125 0.250 0.375 0.500 0.625 0.750 0.875 1.00 

x/c 

Figure 6.17: Circular arc airfoil: Displacement Thickness 

separation point differ by about 5% chord. It is not clear which of the 
methods is in error, but we note that the Cebeci-Smith turbulence model 
is based on equilibrium boundary layer assumptions, and is not valid in or 
near separated zones. 
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^ X/C 

Figure 6.18: Circular arc airfoil: Skin Friction 
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Chapter 7 

CONCLUSIONS 


A finite volume method is applied to the compressible boundary layer 
equations. The novelty of our approach is the use of the integral form 
of the equations without transformation. Instead of incorporating similar- 
ity features into the governing equations via laborious transformations and 
scalings, we use similarity principles explicitly to generate the grid and cal- 
culate a starting solution. Instead of transforming the equations to remove 
the streamwise variation of flow variables in similar flow, we recognize the 
variation and correct for it explicitly when appropriate. 

Important features of the method include: 

• A simple untransformed stream function is used to simplify the gov- 
erning equations. 

• Similarity principles are used explicitly to generate a suitable compu- 
tational grid, to calculate an initial solution, and to construct accurate 
flux formulae. 

• The discrete finite volume equations are linearized and solved via New- 
ton’s method. 

• Reynolds stress terms are calculated with the Cebeci-Smith algebraic 
two equation model. 

• Diverse boundary conditions can be applied easily; the method does 
not distinguish between direct and inverse modes of calculation. 

Results for similar flows and for compressible and turbulent test cases show 
excellent agreement with tabulated solutions and solutions obtained with a 
conventional method. 
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The explicit use of similarity is much simpler than the use of conven- 
tional transformations, but gives equivalent results. Our method retains the 
conceptual simplicity of the integral conservation equations without loosing 
the eifficiency of conventional schemes. Extension to coupled inviscid-viscous 
calculations and to 3-dimensional flows (without the primitive variable for- 
mulation) should prove relatively straightforward. Thus, we recommend the 
method to potential users. 
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Appendix A 

EXACT FLUX 
CALCULATION 


In similar flow , we know the power law behaviour of our variables along 
rj = const. All variables are of the form fc(x) = although <f> may differ 
from variable to variable. We can apply this knowledge to calculate the 
exact fluxes through the streamwise grid lines. 

Straightforward integration of k{x) from X{ to Xi+i gives the exact flux: 



^ + 1 ^ kiXi ), 


l*t 


while our approximation to A.l 


(A.1) 


1**^ — *5 • + ki) • (xj+x Xj), (^-2) 


is accurate to second order in Ax. 

To examine the magnitude of the error define a flux correction factor act 


exact flux __ ) 

approx, flux + ^* ) ’ ( ^t+i ~ ^ 


(A.3) 


For convenience, let ka = and k^ = If = 1, the approximate 
flux, kay is equal to the exact flux, fcje). The larger the deviation of from 
1 the larger the error. Since (A.2) is a linear approximation, ac = 1 only if 
fc(x) is a constant or linear function of x (^ = 0 or 1). 
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Formula A.3 is more instructive if we eliminate x* + Ax*), replace 

ki by axf, and manipulate the resulting expression to get 

_ 2 Xj {1 + Axi/xi)^-^^ - 1 

<f>+lAxi (1 + AX|/xt)^ + 1 

Thus, for a given ^ 0 or 1, ac is a function of Axt/X| only , and, for ^ > 0, 

is bounded by: 

0 < «, < ^ (A.5) 

If Axi/xi is large (and ^ ^ 0, 1) 1 — Oj is large; eis Ax</x,- — *• 0,ac 1- 
We discuss specific applications of ae, to fiuxes calculated for incompressible 
stagnation and flat plate flow, below. 

A.l APPLICATION TO STAGNATION FLOW 

For stagnation flow, simple averaging of the state vectors at i and i + 1 
to get the flux vectors is excellent, because all of the flow quantities vary 
linearly with x along constant q . Since = 1, for both u and r, our 
approximate flux calculations are exeict: ote = 1, identically. 

A.2 APPLICATION TO FLAT PLATE FLOW 

In flat plate flow, u is constant along rj = const., and r ~ x~^!^. Since 
C^e is constant, ^ = 0,«c = 1 and our averaging will yield the exact 

fluxes and uj. However, t\ and Tz will contain some error, since now 
<f> = at = “1/2. Table A.2, below, gives ae and the wall shear stress 
calculated for three values of Ax/x and two different grid densities. r„ 
and Te are the wall shear stresses calculated with ka and ke, respectively. 
For Axi/xi = 10 and Ax,/x< = 1, the truncation error due to equation 
A.2 is large and causes significant error in the predicted value of rw For 
Axi/xi <0.1, (A.2) is so accurate that correction is hardly necessary. 

A cautionary note on the application of equation A.3 is in order. Clearly, 
potential benefits of its application occur only if 

Axi/xi >>~ 0.1 . (A.6) 

However, application of (A.3) for situations where (A.6) holds guarantees 
only that the resulting solution is an accurate approximation to a similar 
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Table A.l: Wall Shear using approximate and exact flux calculations 


Axi/ Xi 

O^c 

^ cells 



10 

0.7120 


0.28122 

0.33328 



■1 

0.28025 

0.33237 

1.0 

0.9706 

MM 

0.32834 

0.33328 



Mm 

0.32744 

0.33237 

0.1 

0.9994 

MM 

0.33319 

0.33328 



mm 

0.33227 

0.33237 




Exact: ; 

= 0.33206 


flow at that station. If the flow is not similar to begin with, one is better off 
using the second order approximation without correction. Another potential 
application of (A.3) is near the leading edge in a flat plate type flow, since 
Xi may be quite small. This use is to be avoided, for it gives the mistaken 
impression that the solution has great physical merit. In fact, the flat plate 
boundary layer approximations are valid only for i2e* >> 1. Solutions 
of the complete Navier-Stokes equations near a leading edge indicate that 
profiles near x = 0 deviate significantly from boundary layer profiles. For 
these reasons, and for simplicity, use of (A.3) for nonsimilar calculations is 
of dubious benefit and is not recommended. 
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Appendix B 

SUTHERLAND’S LAW 


Sutherland models the dependence of viscosity on temperature with 


Mo \To) r + Ti ■ 


(B.l) 


//o is the viscosity at a reference temperature Tq. Ti is a constant: for air 
Ti w 110 K. 

Equation B.l may be written in terms of the enthalpy h = Cp^Ty 


MO 


( ^0 + hi 

\ho) h + hi ^ 


(B.2) 


by eissuming a constant specific heat. 
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Appendix C 


CEBECI- SMITH 
TURBULENCE MODEL 


We use the Cebeci-Smith two equation model to give the eddy viscosity, 
IXi, We outline the relevant equations here, and recommend the thorough 
discussion of the model in [Cebeci 74], 

The model is based on PrandtPs mixing length theory or, equivalently, on 
the Boussinesq eddy-viscosity concept. The Boussinesq formulation models 
the Reynolds stress in the form 

Tt = -/>uV = P^m^- (C.l) 

The Cebeci-Smith model describes the variation of the eddy viscosity in 
the inner and outer regions of turbulent boundary layers by the formulae: 


em.- =L^ 


du 


dy 

€mo=OtUeSl'ftr 


'Ur 


0<y <Vc 
Vc < y < y« , 


(C.2) 


where L is a mixing length, is the kinematic displacement thickness, 'jtr 
is an intermittency factor, and a is a function of the Reynolds number, t/c 
is the normal location where = ^nto- 

In the inner layer, the mixing length L is proportional to the distance 
y from the wall, L = ny. The formula for the inner region is accurate in 
the sublayer and buffer layer if the Van Driest modification for the mixing 
length is used: 

K = K,{l-e-yf% (C.3) 
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where the Karman constant k„ is generally taken to be equal to 0.4, and A 
is the Van Driest damping parameter as modified by Cebeci. 

A takes pressure gradients and/or mass transfer into account: 


A- A+- 



- 1/2 

■ P ■ 

.Pw. 


.Pw . 


nl/2 


where A'*' = 26, and 


N = 

1 - 11.8 


'El]\+ 


• 

.Pe J 

.Pwl 

is a pressure gradient 

P+ = 

\l'eUe 

dUo 


[ur 

dx ’ 


and Ur IS the friction velocity 


(C.4) 


(C.5) 


(C.6) 


(C.7) 


The above formulae used with constant a give accurate results for large 
Reynolds numbers. For small Reynolds numbers agreement with experimen- 
tal results is improved by introducing a variable a. For < 5000, 

ot = f{R$^) as follows: 


a = a<, 


1 + 

1 +3T 


(C.8) 


JTo = 0.55 


= 0.0168 


where 


and 


7T = no [l- exp{-0.2AZz\l^ - 0.298^i)] , 




_ i 


425 


(C.9) 

(C.IO) 

(C.ll) 


The intermittency factor 7 may be user prescribed or experimentally 
defined. 
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C.l NON-DIMENSIONALIZATION 


We non-dimensionalize equations C.2 to C.ll, above, to be consistent 
with the non-dimensionalization introduced in Section 2,2. The form of 
the non-dimensioned equations turns out to be identical to the form of the 
dimensioned equations (a consequence of the dimensional consistency of the 
model and the ‘^natural” non-dimensionalization). We demonstrate this 
below. 

It is convenient to start with the inner model. In non-dimensioned quan- 
tities (simply plugging in relations 2.8), the shear velocity is 




1/2 Fr' 

* w 

■ Pw. 


(C.12) 


As before, primes (’) denote non-dimensioned quantities. Equation C.12 
suggests a non-dimensioned friction velocity 








Now becomes 


P+ = [UooLc„]K [U^]Ul 

p+' = p+ - 


C/o 


[UccU^f [L^ 


dx^ 




dx* 


N is also dimensionless: 



* 



r 1 

2 1 

N = N'^ 

1 - 11.8 



Pe 

p+' 



[p'e\ 


-Pw. 



1/2 


whereeis A becomes 





■ 1 ■ 


fp'l 

N> 

[UooU^l 


-Pw. 



■ 1 ■ 

— 1 
— * 

N' 

kJ 

Ip'wi 


1 1/2 


consistent with 


L'=— = kW 


(C.13) 

(C.14) 

(C.15) 

(C.16) 

(C.17) 

(C.18) 

(C.19) 
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and 


k' = K = Ko(l - e . 

(C.20) 

The coefficient a in equation C.8 is dimensionless; thus, 
equations C.2. These become 

'ffe consider 


(C.21) 

and 

em. = lC^ocLoo}ct[/^S’,;'rtr, 

(C.22) 

so that 


I' ^ \dv-'/dy'\'itr 0<y <ye 

UooLoo \ OtUl6'i*nitr ye<y<Ve 

(C.23) 


As promised, equations C.23 are identical in form to equations C.2. 
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Appendix D 

GLOBAL UNKNOWNS 


This appendix describes the process of including (and solving for) global 
unknowns in the system of equations. Recall that the solution vector is of 
the form 

Auj = dj^i — bj^iA.9^ *“ 

Atpj = ay, 2 - bj^2^9k - Cj,2AI7e 

Ary = ay, 3 - 6y,3A0jfc - Cj^zAUe (151) 

AHj — ay,4 — fey,4A0jt — Cj^^AU^ 

Aqj = ay, 5 — 5y,5A0Jk — Cj^^AU^ 

The influence coefficient vectors a, by and c are known. A9k and AUg are 
global unknowns that must be determined by specifying two auxiliary rela- 
tions in the form: 

AiA9k + BiAUg = i?i . X 

A2A9jt + B2AUg =i?2 ^ ^ ^ 

One of these relations must result from specifying a fifth boundary condition, 
the other can be derived from the definition of yg. The derivation of equation 
D.2 for both cases is discussed below. 

Consider the calculation of the grid edge yg as part of the solution at 
each streamwise station. We wish to impose a condition that ensures that 
yg(x) is always greater than the boundary layer thickness. In most boundary 
layer flows, the kinematic (incompressible) momentum thickness 9k is a good 
measure of the boundary layer thickness [Drela 83]. Thus we define 

y« = nJk , (D-3) 

where r/g is a constant^ chosen such that yg > 6. 

‘Generally = 15, 18 for laminar and turbulent flows, respectively, is adequate. 
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Ok is defined as: 


or, using equation D.3, 


/•»« u 

_ ^ 
~Jo U~e 




1 - 


U 


dy 


drj 


(D.4) 


Equation D.4 implies equation D.3. Thus, we include it as an auxiliary 
relation in the Newton system of equations. 

Newton linearization of eqn. D.4 proceeds as before. We need 


dRese 


dRes0^ 


Reset, + H 5 — ^Awy + ... = 0 . 

ouj^i dui 


(D.5) 


Rese^ = 

}=i 


_ («j+i + ^j) 


1 - 




and 


2 Ue 


2U, 


(»7y+i - Vi) - 1 


“j+l + tij 


Discretizing D.4 gives: 

2Ue 

dResfii^ _ dRese^ _ 
duy dUj+i 

Inserting D.6 - D.7 into D.5 gives 

E(A„,«+A„,)(i-!!i±l±^)^ = E 

i=i * )=i 

which, upon transposing and defining Uati. = -5(uy+i + uy), simplifies to 


2Ui 


Aqy 


(D.6) 


(D.7) 


Uy+i + Uy' 

• Uy+iH-Uy' 

2 

L 2t^e J 


Aqy-Ue 


J-l 

[(Auy+i + Auy)(l/2 - UavelUe) ~ «ave(l “ «o»e/E^e)] Aqy = -Uf 
i=i 

(D.8) 

Equation D.8 is the first of the two additional relations needed to close 
the system of equations. However, the 6uj 's in D.8 are unknown and must 
be related to the global unknowns. We eliminate Auy+i and Auy by using 
relations D.l from the factorized system of equations: 


Auy = ay,i - 6 y,iA^fc - Cy.iAt/^e 
Auy+1 = ay+1,1 - bj+i^iAOk - Cj^i^AUe 


(D.9) 
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Plugging these into Eqn. D.8 and moving known quantities to the right 
hand side gives the required form D.2: 


yZ ” ^ave/Ue)Ar}j} = 


i=i 


j-i 


Ue + ^ “ ^avt/Ut) Wau«(l ^avt/U^] Al]j} 

i=i 


or 


J-l 


^ave/Ue)Ar)j 


i=i 

J-1 


(D.IO) 


*^1 — ^ave/Ue)^Vj 

i=i 


and 


J-i 

~ ^ave/Ue) “ tZavc(l ^ave/^e)] ^Vj 

J=1 

(D.ll) 

The second equation in D.2 is obtained by imposing any one of a number 
of possible boundary conditions. Here we discuss specifying edge velocity, 
wall shear, or displacement thickness; however, other boundary conditions 
can be eeisily implemented. 


1. Specified edge velocity (direct mode): 


Ue - AUe = Uepec 

where Uspec is the specified edge velocity. Thus 


A2 = 0 

B2 = -1 (D12) 

i?2 = Ugpec ^ Ue 


2. Specified wall shear (inverse mode): 

T\ — Afi = Tgpec • 


Using equation D.l 

Ari — ay=i,3 ~ “ Cj=x^z^Ue , 
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gives 


f$pee - ’"1 — -aj=l,3 + + Cj=i^z^Ut 


(D.13) 


so that 

A2 = &j=l,S 
B2 = Cj=lfi 

B2 — Tipee —Ti + aj=i,3 

3. Specified displacement thickness (inverse mode): 

S*-AS* = 6,pec 

Displacement thickness is defined 2is 


rve 

?*= / 

Jo 

rvt 

Jo 

= Ve- 


1 - 

1 - 
1 

PeU, 


JJL 

PtUe 
1 


PeUe dy\ 


dy 


(V'e - V'tu) 


Thus, 


- A5* = 8tptc -yj + 


1 


'^sptc Sf*f * yT 

PjUj 

The left hand side of [D.15] is 


86* 


88* 


rpj (V'lu = 0) 


88* 




via the chain rule, with 


d8*/89k = rtj 




J_ 1 8pj 

TT • 


88*/8Uj= „ , 

pjUj [Uj pj8Uj\ 

88*/8iPj = -1/ipjUj) 


Inserting these relations into [D.15] gives 

- -J^^AUj + \ -Aipj = 8* -yj + 


(D.14) 


(D.15) 


(D.16) 


(D.17) 


PjUj 


PjUj 


PjUj 
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Eliminating At/ij via 


AV'j = oy,2 - bj^2^$k - 

allows us to write this in the usual form (D.2) where 

M = -»?J - bj^2/{pjUj) 

Bi = - (V>J [^+JJ§vi]+ C/.2) 

^2 = ^apec -yj + (V’J - aj^2)/{pjUj) 


(D.18) 
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Appendix E 

NEWTON SYSTEM 


This appendix lists vectors and block matrices defined in the Newton 
system for incompressible laminar flow. In general, A's indicate the differ- 
ence of a quantity taken in the counter-clockwise sense. The exceptions to 
this rule is Arf which is defined as - rfj and Ax which is defined as 
Xt+i - X,-. Quantities with numerical subscripts, say H 2 y are defined as the 
average of H at the adjacent node points. Demonstrative examples of these 
conventions are: 


Au2 — 

Ay2 = Oi^iArij 


U 2 = .5(u,+i,,+i + Ui+i,,) 
«3 = •5(u,-,y+i + «i+l,,+l) 


(E.l) 


The linear equations corresponding to each residual equation in the New- 
ton system 5.4 can be written in block tridiagonal form. For an incompress- 
ible model, this is matrix equation E.2 below. Note that the residuals cor- 
responding to a given cell are shifted by two equations with respect to the 
block structure, due to imposition of the boundary conditions at the wall. 
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J computational cells result in J + 1 unknown changes (5 ’s), and J 
conservation equations. Boundary conditions are implemented in the first, 
second and final equations of the block matrix. The conservation equations 
run from equations 3 through 3 * J - 1 and, due to the boundary conditions, 
are shifted with respect to the block structure. 

Below are the matrices and vectors that make up the Newton system. 
The vector of residuals is: 


i?i = 


0 

0 

T 2 AJ /2 - /xAu2 


(E.3) 




Rj = 


(uiAt/>i + usAi/» 3 + Pi(Ayi + Ays) + (rx - 73) Ai) 

+(u2AV’2 + UiAipi + {P2^V2 + P^Ay^)) 

Axp2 - U2Ay2 
rsAys - /iA«2 

(uiAV'i + usAV’ 3 + Pi(Ayi + Ays) + (ri - rs)Ax) 

+(u2Ai/> 2 + U4AV>4 + (PjAys + P4Ay4)) 

At/>2 - U2Ay2 
0 


The Jacobian block matrices A, B, and C are: 


(E.4) 


(E.5) 


A,= 


1 0 
0 1 
/i 0 


0 

0 

.5Ay2 


(E.6) 


A,= 


Aj = 


B,= 


.5(A^3 + AV> 2 ) «2 - «3 —Ax 
— .5Ay2 1. 0 

/i 0 .5Ay2 

.5(Ai/> 3 + At/> 2 ) U 2 ~ «3 —Ax 

— .5Ay2 1. 0 

1 0 0 

.5(A^i + AV’2) ui - U2 Ax 
— .5Ay2 ~1. 0 

0 0 0 


(E.7) 


(E.8) 


(E.9) 
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(E.IO) 


C,= 


0 0 0 
0 0 0 
-/i 0 .5 Aj/2 


Laminar flow has been assumed for simplicitiy. If in fact a turbulence model 
is implemented it is very important to linearize this eis well. 

The vectors G and H are: 



0 


■(P2 - Pi)Avi 


\P2-P{)Ar,i 

Gi = 

0 

G,= 

-UiArfj 

Gj^ 

-u^Arij 


/zAfji 


T2Ar\j 


0 



■ 0 ■ 


-.5u,+i,j((Ayi -h Ays) -b 2Ay2) 

Hi = 

0 


0 


0 


0 


Hj-i 


+ Ays) + 2Ay2) 

0 

C'j.3,1 


Hj = 


— .5ui+i,j((Ayi + Ays) + 2Ay2) + 
0 


(E.12) 


(E.13) 


(E.14) 


E.l INITIAL CONDITIONS 

The above vectors and block matrices must be modified at the initial sta- 
tion. In the following matrices, it is 2 issumed that a solution is sought at face 
4 (East) and variables are extrapolated from face 4 to faice 2 (West) to obtain 
ordinary difference equations^. Thus, the definitions of the stream function, 
shear stress, and enthalpy flux and their corresponding entries,since they 
contain no terms at face 4, are identical to the marching case. Only the en- 
tries corresponding to the momentum and energy equations are slightly dif- 
ferent; and, for this incompressible case, we need give only the first row/entry 
in each of the matrices/vectors. 

Some new definitions are (to implement similarity principles): 

Xu = (xf/xi+i)"*, xt = (x,/x,+i)“‘ 

‘Alternatively, the solution could be sought at face 2 and variables extrapolated from 
face 2 to face 4. 
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, = {Xi/Xi+lY\ Xy = 

®ul3 = 1 + ®u> ®tl3 = 1 + It 
I«13 ~ 1 ~ ®yl3 = 1 ~ ly . 

With these definitions, the following identities are true in similar flow: 

A^l = -xl>ijXgi3, A^3 = ipi,j+lXsl3 
A^2 = -AV’4i«, etc. for Ay 

and 

Ui = .5u,-,yX„i3, U3 = .5 u,-,,+iXu13 

U2 = .5u4Xu, etc. for t, H, q ,tp . 


Rj,i — 


(uiAV»i + uzAips + Pi(Ayi + Ay3) + (ri - r3)Ax) 

+(u 2AV»2 + «4AV’4 + (p2Ay2 + PiAyt)) 


(E.15) 


Bj,rowl = [•5(Zu13AV’1 + AV>2®u + AV»4) -«ll*13 “ «2l» + «4 .SXtlsAxj 

(E.16) 

^jyrowl — \ . 5 {ZulS^ifZ + Atp 2 ^u + ^^4) t*3X,i3 + U2X, “ «4 .Sx^isAx] (E.17) 
Gj^i = (Pi(Ayi + At/3) + + P4Ay4)/0t j (E.18) 

_ r + xj) * (Ayi + Ays) 

[ +2(xjAy2 + Ay4)) 

Clearly, changes in the discrete equations to allow calculation of an ini- 
tial solution are minor and can easily be included 21s a special case of the 
marching solution. The omitted matrices and vectors (at j = 1, j = J — 1, 
and i = J) can be constructed by simple analogy with the marching solution 
setup above. 


(E.19) 
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Appendix F 

PROGRAM LISTING 


A program listing of the FORTRAN computer code is available upon 
request. 
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